!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine write_vtk_unstructured(datavtk_name, var_name, var, x, y, nx, ny, n_blocks_all, n_var, var_dim, format, header)
!==================================================================================================================================|

   implicit none

   integer(4), parameter :: n_dim = 3
   integer(4), parameter :: n_point_cell = 4
   integer(4), parameter :: file_unit = 21

   integer(4), intent(in) :: nx, ny, n_blocks_all, n_var, format

   character(*), intent(in) :: var_name(n_var), datavtk_name
   character(256), intent(in) :: header

   integer(4), intent(in) :: var_dim(n_var)
   real(4), intent(in) :: x(nx, n_blocks_all), y(ny, n_blocks_all)
   real(4), intent(in) :: var(sum(var_dim), nx, ny, n_blocks_all)

   character(12) :: int_num1, int_num2
   integer(4) :: n_point, n_cell
   integer(4) :: i_var, i, j, index_var, i_point, i_block, pos
   integer(4) :: nxp1, nxny
   integer(4), allocatable :: cell_type(:), connectivity(:, :), conn(:, :)
   real(4), allocatable :: point(:, :)

!----------------------------------------------------------------------------------------------------------------------------------|
!  prepare
!----------------------------------------------------------------------------------------------------------------------------------|
   n_point = nx * ny * n_blocks_all
   n_cell = (nx - 1) * (ny - 1) * n_blocks_all
   nxp1 = nx + 1
   nxny = nx * ny

   allocate (cell_type(n_cell))
   allocate (connectivity(n_point_cell, n_cell))
   allocate (conn(n_point_cell + 1, n_cell))
   allocate (point(n_dim, n_point))
   cell_type = 9

   i_point = 1
   do i_block = 1, n_blocks_all
   do j = 1, ny
   do i = 1, nx
      point(1, i_point) = x(i, i_block)
      point(2, i_point) = y(j, i_block)
      point(3, i_point) = 0.0
      i_point = i_point + 1
   enddo
   enddo
   enddo

   i_point = 1
   do i_block = 1, n_blocks_all
   do j = 1, ny - 1
   do i = 1, nx - 1
      pos = (i - 1) * 1 + (j - 1) * nx + (i_block - 1) * nxny
      connectivity(1, i_point) = pos
      connectivity(2, i_point) = pos + 1
      connectivity(3, i_point) = pos + nxp1
      connectivity(4, i_point) = pos + nx
      i_point = i_point + 1
   enddo
   enddo
   enddo

   conn(1, :) = n_point_cell
   conn(2 : n_point_cell + 1, :) = connectivity

!----------------------------------------------------------------------------------------------------------------------------------|
!  binary data
!----------------------------------------------------------------------------------------------------------------------------------|
   if (format .eq. 0) then
      open(file_unit, file = trim(datavtk_name), form = 'binary', convert = 'big_endian')
!----------------------------------------------------------------------------------------------------------------------------------|
!     write the head data
!----------------------------------------------------------------------------------------------------------------------------------|
      write(file_unit) '# vtk DataFile Version 2.0' // char(10)
      write(file_unit) trim(header) // char(10)
      write(file_unit) 'BINARY' // char(10)
      write(file_unit) 'DATASET UNSTRUCTURED_GRID' // char(10)

!----------------------------------------------------------------------------------------------------------------------------------|
!     write the point information
!----------------------------------------------------------------------------------------------------------------------------------|
      write(int_num1, "(i12)") n_point
      write(file_unit) 'POINTS '// int_num1 // ' float' // char(10)
      write(file_unit) point

!----------------------------------------------------------------------------------------------------------------------------------|
!     write the cell information
!----------------------------------------------------------------------------------------------------------------------------------|
      write(int_num1, "(i12)") n_cell
      write(int_num2, "(i12)") (n_point_cell + 1) * n_cell
      write(file_unit) 'CELLS ' // int_num1 // ' ' // int_num2 // char(10)
      write(file_unit) conn
      write(file_unit) 'CELL_TYPES ' // int_num1 // char(10)
      write(file_unit) cell_type
!----------------------------------------------------------------------------------------------------------------------------------|
!     write the point data
!----------------------------------------------------------------------------------------------------------------------------------|
      write(int_num1, "(i12)") n_point
      write(file_unit) 'POINT_DATA ' // int_num1 // char(10)

      index_var = 1
      do i_var = 1, n_var
!----------------------------------------------------------------------------------------------------------------------------------|
!        write the scalar point data
!----------------------------------------------------------------------------------------------------------------------------------|
         if (var_dim(i_var) .eq. 1) then
            write(file_unit) 'SCALARS ' // trim(var_name(i_var)) // ' float' // char(10)
            write(file_unit) 'LOOKUP_TABLE default' // char(10)
            write(file_unit) var(index_var, :, :, :)
            index_var = index_var + var_dim(i_var)
!----------------------------------------------------------------------------------------------------------------------------------|
!        write the vector point data
!----------------------------------------------------------------------------------------------------------------------------------|
         elseif (var_dim(i_var) .eq. 3) then
            write(file_unit) 'VECTORS ' // trim(var_name(i_var)) // ' float' // char(10)
            write(file_unit) var(index_var : index_var + 2, :, :, :)
            index_var = index_var + var_dim(i_var)
         endif
      enddo
      close (file_unit)

!----------------------------------------------------------------------------------------------------------------------------------|
!  ASCII data
!----------------------------------------------------------------------------------------------------------------------------------|
   else
      open(file_unit, file = trim(datavtk_name), form = 'formatted')
!----------------------------------------------------------------------------------------------------------------------------------|
!     writethe head data
!----------------------------------------------------------------------------------------------------------------------------------|
      write(file_unit, '(a)') '# vtk DataFile Version 2.0'
      write(file_unit, '(a)') trim(header)
      write(file_unit, '(a)') 'ASCII'
      write(file_unit, '(a)') 'DATASET UNSTRUCTURED_GRID'

!----------------------------------------------------------------------------------------------------------------------------------|
!     write the point information
!----------------------------------------------------------------------------------------------------------------------------------|
      write(file_unit, '(a, 1x, i, 1x, a)') 'POINTS', n_point, 'float'
       write(file_unit, '(6es)') point

!----------------------------------------------------------------------------------------------------------------------------------|
!     write the cell information
!----------------------------------------------------------------------------------------------------------------------------------|
      write(file_unit, '(a, 1x, i, 1x, i)') 'CELLS', n_cell, (n_point_cell + 1) * n_cell
      write(file_unit, '(9i)') conn
      write(file_unit, '(a, 1x, i)') 'CELL_TYPES', n_cell
      write(file_unit, '(12i)') cell_type

!----------------------------------------------------------------------------------------------------------------------------------|
!     write the point data
!----------------------------------------------------------------------------------------------------------------------------------|
      write(file_unit, '(a, i)') 'POINT_DATA', n_point

      index_var = 1
      do i_var = 1, n_var
!----------------------------------------------------------------------------------------------------------------------------------|
!        write the scalar point data
!----------------------------------------------------------------------------------------------------------------------------------|
         if (var_dim(i_var) .eq. 1) then
            write(file_unit, '(a, 1x, a, 1x, a)') 'SCALARS', trim(var_name(i_var)), 'float'
            write(file_unit, '(a)') 'LOOKUP_TABLE default'
            write(file_unit, '(6es)') var(index_var, :, :, :)
            index_var = index_var + var_dim(i_var)
!----------------------------------------------------------------------------------------------------------------------------------|
!        write the vector point data
!----------------------------------------------------------------------------------------------------------------------------------|
         elseif (var_dim(i_var) .eq. 3) then
            write(file_unit, '(a, 1x, a, 1x, a)') 'VECTORS', trim(var_name(i_var)), 'float'
            write(file_unit, '(6es)') var(index_var : index_var + 2, :, :, :)
            index_var = index_var + var_dim(i_var)
         endif
      enddo
      close (file_unit)
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  free memory
!----------------------------------------------------------------------------------------------------------------------------------|
   deallocate (cell_type, connectivity, conn, point)

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine write_vtk_unstructured
